Table of Contents

New LaTeX commands declared in this cell. Expand to view. $$ \newcommand\abs[1]{\left|#1\right|} \newcommand\dx[2]{\frac{\partial#1}{\partial#2}} \newcommand\ddx[2]{\frac{\partial^2#1}{\partial#2^2}} \newcommand\divv[1]{\nabla\cdot\vect{#1}} \newcommand\curl[1]{\nabla\times\vect{#1}} \newcommand\grad[1]{\nabla\vect{#1}} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\braket}[1]{\langle#1\rangle} \newcommand{\aint}[0]{\int_{-\infty}^{\infty}} \newcommand{\vect}[1]{\mathbf{#1}} \newcommand{\order}[1]{\mathcal{O}{#1}} $$

Fourier Series

See discussions in:

Fourier Transform

Basics

The following discussion assumes evenly spaced steps in time $\Delta t$. The Fourier series of a function is its decomposition into individual sines and cosines:

$$ y(t) = \sum_{j=1}^n y_j \sin(2 \pi f_j t + \phi_j) \\ \text{With } _j = \frac{\omega_j}{2\pi} $$

.

If we take this to the continuum limit, summing over all frequencies,

$$ y(t) = \int_\infty^\infty Y(f) e^{-2\pi ift} df = \frac{1}{2\pi} \int_\infty^\infty Y(\frac{\omega}{2\pi}) ^{-i\omega t} d\omega \\ Y(f) = \int_\infty^\infty e^{2\pi i f t} dt$$

The two expressions are equivalent, though the $\omega$ version is more commonly seen in the context of quantum mechanics. The forms above are commonly called the forward Fourier transforms, though this has no deep meaning. The idea of going 'forward' likely results from our far more pedestrian experience of time as a dimension than frequency; we transform 'forward' to a new domain, and then inverse transform to return to 'normal'. This reasoning holds for the other common Fourier transforms as well; consider that position space is far more familiar to us than momentum space!

Another tool that will be useful is the Dirac delta function, the result of transforming forward, then back, and expecting your original function to return (a reasonable hope):

$$ \int_\infty^\infty e^{i\omega(t - t')} d\omega = 2\pi \delta(t-t')\\ $$

This function has the property of being zero everywhere except where t=t'. It has the following property when integrated with a function:

$$ \int y(t') \delta(t-t') dt' = y(t) $$

It thus has the inverse units of your integration variable.

Discrete FT

Consider discretizing the above definition of the Fourier transform:

$$ \text{For }t = m \Delta t \text{, and} f_n = n/(N\Delta t): \\ y_m = \sum_{n=0}^{N-1} Y_n e^{-2\pi imn/N} \qquad Y_n =\frac{1}{N} \sum_{m=0} y_m e^{2 \pi i m n / N} $$

We can recreate the discrete version of a Dirac delta function with these definitions; the Kronecker delta function:

$$ \sum_{n=0}^{N-1}e^{2 \pi i n(m - m')/N} = N \delta_{m,m'} $$

Note that $\Delta t$ has disappeared from our discrete FT. There are subtleties involving the amount of information contained in the DFT between real and complex cases.

The DFT as implented below loops through the time and frequency indices, m and n, to compute matrix elements $$\bra{n} M \ket{m}$$, where $M$ is the DFT Matrix. Since this matrix has $N\times N$ elements, the computation is of $\mathcal{O}(N^2)$, which will quickly become cumbersome. The Fast Fourier Transform is an algorithm much more efficient at computing the FT, which will be discussed in the next section.

Nyquist Frequency

The Nyquist critical frequency is defined as $$ f_c = \frac{1}{2\Delta t} $$ The simplest case to consider is a sine wave: to properly sample a sine wave, you must take at least 2 points per cycle. This means that you need to sample at a minimum of twice the maximum frequency of your data to preserve information below that frequency. In a real device, such as an experimental setup limited by the response of particular electrical cables, there will probably be natural frequency cutoffs. You need to sample at twice the max frequency you want. For example, human hearing cuts off at around 20kHz. To be able to properly sample audio recordings for future listening (on headphones, for instance), you need to sample at twice this frequency. It's no accident that audio CDs sample at 44.1kHz!

If the datatype is not set to complex, you will be discarding the phase of your DFT; $$ e^{i\theta} = \cos{\theta} + i \sin{\theta} $$ With only half of the information, you will introduce phase shifts into the Fourier transform. This can be easily tested by changing the functions above to read dtype=float.

FFT

The Fast Fourier Transform (FFT) is a ubiquitous algorithm for computing the Fourier transform of a dataset in $\mathcal{O}(N \log N) $ time. Compare the difference in the number of operations needed for the two cases.

$$ N = 8 = 2^3 \\ Y_n = \sum_{m=0}^{N-1} M_{n,m}Y_m \\ M_{nm} = e^{2 \pi i n m / N} \\ ex: Y_0 = \sum_{m=0}^{7} e^{2\pi i (0) m / N}y_m = y_0 + y_1 +...\\ \text{Split into even (e) & odd (o) m:} \\ Y_n = Y^e_n + w^n Y_n^o \qquad w = e^{2 \pi i / N} $$
$$ = \sum_{m'=0}^3 Y_{2m'} w^{2m'n} + w^n \sum_{m'=0}^3 Y_{2m'+1} w^{2m'n} \\ Y_n = Y^{e}_{n} + w^n y_n^o $$

For each level, we can recursively split the evens and odds again and get twice as many terms:

$$ Y^e_n = Y^{ee}_n + w^{2n} Y^{eo}_n \\ Y^o_n = Y^{oe}_n + w^{2n} Y^{oo}_n \\ $$

split until each $Y^{eoeoeoeoeoooeooeo...eoeooooeooooeeeeeoeoee}_n$, etc, term is a single element. For the n=8 example, this means three iterations of dividing the input array.

$$ Y^{ee}_n = Y^{eee}_{(0)} + w^{4n} Y^{eeo}_{(4)} \\ Y^{eo}_n = Y^{eoe}_{(2)} + w^{4n} Y^{eoo}_{(6)} \\ Y^{ee}_n = Y^{ooe}_{(3)} + w^{4n} Y^{ooo}_{(7)} \\ Y^{ee}_n = Y^{oee}_{(1)} + w^{4n} Y^{oeo}_{(5)} \\ $$

Recall what these Y terms are. Each is a broken down component of the initial array $Y_n$, $n \in [0,8).$ A string like eeo, eoe, etc, represents a reversed binary representation of the index to $Y_n$. $e\rightarrow 0, o \rightarrow 1$. If you flip the text string, you recover the binary representation of the index. An example is essential here.

$$ Y^{eoo}: eoo \rightarrow 011; \text{reversed, } 110 \rightarrow Y_6 \\ Y^{eoo} = Y_6 $$

So we can bit-reverse the index of the original array to find the new location of that element after the FFT algorithm begins working.

$$ Y_m = \sum_{n=0}^{N-1} y_n e^{-2\pi i m n / N} \\ $$ $$\order{N \log_2{N}}$$

Bit reversal

The point of bit reversing is a bookkeeping shortcut. Rather than splitting down the levels, sorting even/odd elements until they come out as a new array

$$[y_0, y_4, y_2, y_6, y_1, y_5, y_3, y_7]$$

There is a simple pattern here that is not obvious at first glance. We can recognize that, if we convert the indices to binary, with the number of bits corresponding to $ln_2(len(y))$, we get

$$[y_{000}, y_{100}, y_{010}, y_{110}, y_{001}, y_{101}, y_{011}, y_{111}]$$

And, if we reverse the order of each of these bits, the final result is

$$[y_{000}, y_{001}, y_{010}, y_{011}, y_{100}, y_{101}, y_{110}, y_{111}]$$

a slick trick that can save us time by instantly swapping our input array elements to the place that sorting would have put them already.

Recursive FFT

Below is a sketch of the recursive fft algorithm. Note that this does not utilize the bit reversal technique; deriving this method led to the recognition of bit reversal as a bookkeeping trick, which is implemented in Iterative FFT below.

def _recursive_fft(data):
    """recursively compute the fft of the input data."""
    N = int(len(data))
    if N == 1:
        return data # The FT of a single element is itelf
    w_N = np.exp(-2j * np.pi / N)
    w = 1.0 + 0.0j

    data_even = data[0::2] # starts from 0 and moves two elements until N-2
    data_odd = data[1::2] # starts from 1 and moves two elements until N-1

    y_e = recursive_fft(data_even)
    y_o = recursive_fft(data_odd)

    y = np.zeros([N], dtype=complex)
    for k in range(0, N//2):
        t = w * y_o[k] # butterfly operation
        y[k] = y_e[k] + t
        y[k+N//2] = y_e[k] - t
        w = w * w_N
        del t
    return y

Iterative FFT

The above recursive FFT is messy (I'm not sure where the extra frequency in the fft is coming from), and padding zeros is difficult without significantly more bookkeeping than I've done above. We will move from the bottom up in the tree describing the FFT, following the bit-reversal scheme.

def iterative_fft(cls, data):
    data = data.astype(complex)
    assert (np.log2(len(data)) % 1) == 0
    A = cls.reverse_array(data)
    N = int(len(A))
    S = int(np.log2(N))
    for s in range(1, S+1):
        m = int(2**s) # length of split array at each level
        w_m = np.exp(-2j * np.pi / m)
        for k in range(0, N, m): # by m
            w = 1.0 + 0.0j
            m_2 = (m // 2)
            for j in range(0, m_2):
                t = w * A[k + j + m_2]
                u = A[k + j]
                A[k + j] = u + t
                A[k + j + m_2] = u - t
                w = w * w_m
        return A

Bit reversal

def bit_reverse(num, bits):
    num = int(num)
    return int('{:0{s}b}'.format(num, s=bits)[::-1], 2)

Comparative Efficiency

For this simple problem, poorly implemented, recursive_fft() runs 38 times faster than my $\mathcal{O}(N^2)$ DFT routine, while numpy's far superior np.fft.fft runs 2292 times faster than recursive_fft().

Sample interval and N

The interval over which to sample and the number of samples are the parameters most closely in our control. It was briefly mentioned above that discarding the imaginary part of the fft destroys the phase information, but this may have been obfuscated by the development of the crude fft. Let's see how a phase shift in a pure $\sin{\omega t}$ signal affects the fourier transform.

Good frequencies are $$ f = \frac{n}{T} \bigg \rvert_{-N/2}^{N/2 - 1} \\ \omega = \frac{2\pi n}{T} \bigg \rvert_{-N/2}^{N/2 - 1} \\ k = \frac{2\pi n}{L} \bigg \rvert_{-N/2}^{N/2 - 1} \\ $$

From the first row of plots, you should be able to convince yourself that there should only be an imaginary piece to the frequency spectrum. From euler's equation, $$e^{i\omega t} = \cos{\omega t} + i \sin{\omega t},$$ the FT depends on a sine and a cosine contribution, but our initial signal was $\sin{\omega t}$. The cosine amplitude should then be zero.

In the second pair of plots, a phase shift was applied. The phase shift can be deconstructed into a relative weighing of a sine and a cosine piece, with relative amplitudes. This is what we see; there is now a real and an imaginary amplitude after the FFT, both at the same frequency but with different amplitude. This was the cause of the loss of phase and amplitude earlier; without specifying that the datatype was complex, Python dumped the complex part of the arrays to cast to real floats, and we lost our phase and amplitude information.

In the plot below, frequencies incommensurate with the sampling frequency were picked (and deliberately chosen to be primes, though for no specific reason). Even though only two frequencies (5Hz and 7Hz) are present in the signal data, there are nonzero amplitudes at multiple frequencies, with the maximum amplitude being where you'd expect them; i.e., at the given frequencies in the summed sine waveforms.

The frequency spectrum has twice as many peaks as you'd expect; the components $Y_{N/2-n} = Y^*_{N/2+n}$. Since $ N/2 \rightarrow F_{\text{nyquist}} $, the duplicated frequencies will be rotated about the Nyquist frequency.

Aliasing

The above discussion only used frequencies below the Nyquist frequency. If we sample above this, we run into problems.

Let's look at a signal of 10Hz, with the sampling chosen so that we have a nyquist frequency of 6Hz.

The doubled fft frequency spectrum contains what looks like a 10Hz signal, but if we look at the domain [0,f_nyq] it would appear that there is a signal at 2Hz. If we were handed this spectrum blind, we would have no way of knowing if the 'real' frequency of the input signal was the 2Hz or 10Hz peak. This is aliasing.

The full domain is from 0 to the sampling frequency, and the doubled signals are obtained by folding the plot from the center. The difference between the input frequency and the nyquist frequency is $$ f_0 - f_{nq} = 10Hz - 6Hz = 4Hz $$

If we move that distance from $f_{nq}$ again, we'll find our aliased frequency.

$$ f_{nq} - (f_0 - f_{nq}) = 2f_{nq} - f_0 = 2Hz $$

This is what we see in the plot above. Let's look at the time domain.

Wagon Wheel effect on Wikipedia

The fourier transformed data comes back with an aparrent frequency far lower than the input frequency, since we sampled lower than half the data's frequency. In handling real data one would want to either have a way of estimating the frequency before sampling, or taking as high a sample rate as feasible, before analysing the data with signal processing methods.

Derivatives

We can use the Fourier transform to take derivatives of analytic, periodic functions.

$$ f(x) = e^{\cos{2\pi x/L}} \\ f'(x) = -\frac{2\pi}{L}f(x) \sin\frac{2\pi x}{L} \\ f''(x) = \frac{2\pi}{L}\bigg[- f'(x)\sin\frac{2\pi x}{L} - f(x)\frac{2\pi}{L}\cos{\frac{2\pi x}{L}} \bigg]\\ = \frac{4\pi^2}{L^2} f(x) \bigg[\sin^2{\frac{2\pi x}{L}} - \cos{\frac{2\pi x}{L}} \bigg] $$

Errors

Errors come in two flavors from physics lingo - UV and IR.

Ultraviolet (UV)

A wave cannot fit in a box with a frequency that will touch less than two points. I.E., the smallest wavelength that can fit in your box is 2*dx. Imagine an 8 point lattice with a standing wave of this wavelength:

$$ \text{(1) } \uparrow \downarrow \uparrow \downarrow \uparrow \downarrow \uparrow \downarrow \\ \text{(2) } \uparrow \uparrow \uparrow \uparrow \downarrow \downarrow \downarrow \downarrow $$

Case (1) has the smallest wavelength possible. Anything less than this and we can't represent it accurately; this will cause aliasing, and the resultant aliased wavelength will be larger than you want (lower energy). See Numerical Recipes for a discussion of aliasing as energy bleeding out of the edge of the PSD power spectrum.

These are called ultraviolet errors because the low wavelengths correspond to high energy waves (particles); thus the errors are in the ultraviolet regime on the electromagnetic scale.

Infrared (IR)

Likewise, a wave with too large a wavelength will begin to lose features. If the wavelength is larger than case (2) above, you may end up with IR errors. We could fit a half wavelength into a box, but our periodic boundary conditions will make a discontinuity appear (consider what happened when there was just one extra zero in a Fourier transformed data set from improperly picking boundary conditions!)

Gaussian Example

Consider a Gaussian in a centered box:

$$ f(x) = e^{-x^2 / 2 \sigma^2}, \qquad f'(x) = \frac{-x}{\sigma^2}f(x), \qquad f(k) = \sqrt{2\pi} \sigma e^{-k^2 \sigma^2/2} $$

We want the function to achieve machine precision. In this case, that means approaching $\epsilon = 2^{-52} \approx 2 \times 10^{-16}$.

This estimation scheme for optimal values of N and L is great if we have an analytic form to play with, but what if the function is too complicated for this? We have a few options:

Power Spectrum (incomplete)

Obtaining the power spectrum via fft requires using the convolution of two functions.

Convolution

Correlation / Autocorrelation

We will need to integrate a function to obtain the correlation between it and itself, lagged. The correlation is defined as $$ Corr[y](\tau) = \int_{\infty}^{\infty} y(t)*y(t + \tau)d\tau $$

Correlation: $$ Corr[A, B][n] = \sum_{m = 0}^{m= N-1} A[m]B[n-m] $$

Parseval's Theorem

Parseval's Theorem relates the total energy in the time domain to the total energy in the frequency domain. The Fourier transform takes you from one basis to another, but the physical quantities associated with the system cannot change. Consider this in the continuous and discrete cases:

$$ \int dt e^{-i\omega t} y(t) = Y(\omega) \qquad \frac{1}{2\pi} \int d\omega e^{i \omega t} Y(\omega) = y(t) \\ \sum_{n=0}^{N-1} e^{-i \omega n / N} y_n = Y_m \qquad \frac{1}{N} \sum_{m=0}^{N-1} e^{i \omega m / N} Y_m = y_n \ $$

Parseval's Theorem states that:

$$ \boxed{ \int dt \abs{y(t)}^2 = \frac{1}{2\pi}\int d\omega \abs{Y({\omega})}^2 \\ \sum_{n=0}^{N-1} \abs{y_n}^2 = \frac{1}{N} \sum_{n=0}^{N-1} \abs{y_m}^2 } $$

Numerical Issues

$$ \nabla^2(ab) = \nabla (b \nabla a + a \nabla b) \\ = \nabla b \cdot \nabla a + b \nabla ^2 a + \nabla a \cdot \nabla b + a \nabla ^2 b \\ \rightarrow \nabla a \cdot \nabla b = \frac{\nabla ^2 (ab) - a \nabla b - b \nabla a}{2} \\ $$

Why does this fail? The set of frequencies accessible to a and b are obtained with np.fftfreq(), but a*b has a different set of frequencies. Borrowing from a discussion by MMF, the dft has frequencies

$$ a_k = \sum_k e^{-ikx}a_x, \qquad k_n = \frac{2\pi n}{L} \in \bigg [ \frac{-N}{2} \leq n \leq \frac{N}{2} \bigg ) $$

If we multiply our two functions $a_x$ and $b_x$, the spectrum now looks like

$$ a_x b_x = \frac{1}{N^2} \sum_{k_a, k_b} e^{i(k_a+k_b)x}A_{k_1}B_{k_2} = \frac{1}{N}\sum_p e^{ipx}C_p $$

In other words, while the analytic derivative satisfies the product rule as used above, the FFT derivative will not be able to 'see' the frequencies the new function have outside the old frequency basis.

Vector Basis

$$\begin{gather} \bra{x}\ket{g} = g(x) \\ \bra{x}\ket{f} = f(x) \\ \bra{f}\ket{x}\bra{x}\ket{g} = f^\ast(x)g(x)\\ \int_{-\infty}^{\infty}dx \bra{f}\ket{x}\bra{x}\ket{g} = \bra{f}\ket{g} \\ = \int_{-\infty}^{\infty}\frac{dk}{2\pi} f^\ast(k)g(k) \\ \bra{k}\ket{x} = e^{-ikx} \\ \rightarrow \int\frac{dk}{2\pi}e^{ikx}f(k) = \int\frac{dk}{2\pi}\bra{x}\ket{k}\ket{k}\bra{f} \\ = \bra{x}\ket{f} = f(x) \\ \end{gather} $$

Example: Filter a signal with an FFT

Links: